# ------------------------------------------------------------------------------
# Tamaki & Jung 2025
# The Non-linearity of Populist Attitudes and Political Ideology
# Replication Code for PSRM
# Supplementary Materials 
# ------------------------------------------------------------------------------
# Log file 
sink("Tamaki_and_Jung_replication_Appendix_log.txt", split = TRUE)
setwd("C:/Users/yujin/Dropbox/Cowork-Populist Attitudes/PSRM - R and R/Final Code")
# PREAMBLE ---------------------------------------------------------------------
# Clean Environment
rm(list = ls())
gc()

# Load Required Packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, mgcv, sandwich, lavaan)

# Load Data
load("df_final(11.24).RData")

# ------------------------------------------------------------------------------
# Data Preparation
# ------------------------------------------------------------------------------
df2 <- df %>%
  select(pop_add_m, pop_add, ideo, ext2,
         age, gen, edu, empl,
         corrupt_ind, econ_ind, pol_int,
         socecon_status, occ_status,
         country_name) %>%
  mutate(
    country_name = factor(country_name),
    socecon_status = case_when(
      socecon_status == 1 ~ "White Collar",
      socecon_status == 2 ~ "Worker",
      socecon_status == 3 ~ "Farmer",
      socecon_status == 4 ~ "Self-employed",
      TRUE ~ NA_character_
    ),
    socecon_status = relevel(factor(socecon_status), ref = "White Collar")
  )

# ------------------------------------------------------------------------------
# Appendix A – Survey Demographics 
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Appendix B – Measurement Invariance Tests (MGCFA)
# ------------------------------------------------------------------------------
# 1. Prepare Data for CFA
pop_var <- df %>%
  select(m1, ae1, ae3, ae4, pc1, country_name) %>%
  na.omit() %>%
  filter(!country_name %in% c(
    "Greece_2015", "Hong Kong", "Ireland",
    "Republic of Korea", "Taiwan_2016"
  ))

# 2. Specify Measurement Model (excluding AE2 due to reverse wording)
pop_cfa <- "pop_n =~ m1 + ae1 + ae3 + ae4 + pc1"

# 3. Configural Invariance -----------------------------------------------------
fit_mgcfa_conf <- cfa(
  model = pop_cfa,
  data = pop_var,
  estimator = "mlr",
  missing = "fiml",
  group = "country_name"
)

# Save configural summary to .txt
summary_conf <- capture.output(summary(fit_mgcfa_conf,
                                       fit.measures = TRUE,
                                       standardized = TRUE))

# Save configural summary to .txt
if (!dir.exists("Data")) dir.create("Data")
writeLines(summary_conf, con = "Data/(01) MGCFA_Config_Summary.txt")

# 4. Metric Invariance ---------------------------------------------------------
fit_mgcfa_metric <- cfa(
  model = pop_cfa,
  data = pop_var,
  estimator = "mlr",
  missing = "fiml",
  group = "country_name",
  group.equal = "loadings"
)

# Save metric summary to .txt
summary_metric <- capture.output(summary(fit_mgcfa_metric,
                                         fit.measures = TRUE,
                                         standardized = TRUE))

if (!dir.exists("Data")) dir.create("Data")
writeLines(summary_metric, con = "Data/(02) MGCFA_Metric_Summary.txt")

# 5. Compare Configural vs Metric Models ---------------------------------------
anova_mgcfa <- anova(fit_mgcfa_conf, fit_mgcfa_metric)

# ------------------------------------------------------------------------------
# Appendix C.1 – Stepwise Models: Populist Attitudes ~ Ideology
# ------------------------------------------------------------------------------
# Model 1: Main Effect Only (Smooth term for ideology)
gam_m1 <- gam(
  pop_add_m ~ s(ideo) +
    factor(country_name),
  method = "REML",
  data = df2
)

# Model 1.2: + Basic Demographic Controls
gam_m1.2 <- gam(
  pop_add_m ~ s(ideo) +
    age + gen + edu + empl +
    factor(country_name),
  method = "REML",
  data = df2
)

# Model 1.3: + Full Controls
gam_m1.3 <- gam(
  pop_add_m ~ s(ideo) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df2
)

# Model Comparison: AIC and BIC
AIC(gam_m1, gam_m1.2, gam_m1.3)
BIC(gam_m1, gam_m1.2, gam_m1.3)

# Summaries
summary(gam_m1)
summary(gam_m1.2)

# Table Output: Ideology Stepwise Models
stargazer::stargazer(
  gam_m1, gam_m1.2, gam_m1.3,
  type = "text",
  out = "Table_C1.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001",
  notes.append = FALSE
)

# ------------------------------------------------------------------------------
# Appendix C.2 – Stepwise Models: Populist Attitudes ~ Ideological Extremism
# ------------------------------------------------------------------------------
# Model 2: Main Effect Only (Smooth term for extremism)
gam_m2 <- gam(
  pop_add_m ~ s(ext2) +
    factor(country_name),
  method = "REML",
  data = df2
)

# Model 2.2: + Basic Demographic Controls
gam_m2.2 <- gam(
  pop_add_m ~ s(ext2) +
    age + gen + edu + empl +
    factor(country_name),
  method = "REML",
  data = df2
)

# Model 2.3: + Full Controls
gam_m2.3 <- gam(
  pop_add_m ~ s(ext2) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df2
)

# Model Comparison: AIC and BIC
AIC(gam_m2, gam_m2.2, gam_m2.3)
BIC(gam_m2, gam_m2.2, gam_m2.3)

# Summaries
summary(gam_m2)
summary(gam_m2.2)

# Table Output: Extremism Stepwise Models
stargazer::stargazer(
  gam_m2, gam_m2.2, gam_m2.3,
  type = "text",
  out = "Table_C2.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001",
  notes.append = FALSE
)

# ------------------------------------------------------------------------------
# Appendix D: Political Interest
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# (1) Create Political Interest Category Variable
df2 <- df2 %>%
  mutate(
    pol_int_cat = case_when(
      pol_int == 1 ~ 1,                    # Very interested
      pol_int == 2 ~ 2,                    # Interested
      pol_int %in% c(3, 4) ~ 3,            # Not interested
      TRUE ~ NA_real_
    ),
    pol_int_cat = as.factor(pol_int_cat)
  )
# ------------------------------------------------------------------------------
# (2) – Model 1: Ideology + Political Interest as Control
pi.gam_m1 <- gam(
  pop_add_m ~ s(ideo) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(pol_int_cat) +
    factor(country_name),
  method = "REML",
  data = df2
)

summary(pi.gam_m1)

# Diagnostics
gam.check(pi.gam_m1)               # Smooth term diagnostic (p > 0.05)
concurvity(pi.gam_m1, full = TRUE) # Max concurvity for ideology ≈ 0.09 (acceptable)
# ------------------------------------------------------------------------------
# (3) – Model 2: Extremism + Political Interest as Control
pi.gam_m2 <- gam(
  pop_add_m ~ s(ext2) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(pol_int_cat) +
    factor(country_name),
  method = "REML",
  data = df2
)

summary(pi.gam_m2)

# Diagnostics
gam.check(pi.gam_m2)               # Smooth term diagnostic (p > 0.05)
concurvity(pi.gam_m2, full = TRUE) # Max concurvity for extremism ≈ 0.09 (acceptable)
# ------------------------------------------------------------------------------
# (4) – Table Output for Models with Political Interest
stargazer::stargazer(
  pi.gam_m1, pi.gam_m2,
  type = "text",
  out = "Table_D.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "pol_int_cat", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001",
  notes.append = FALSE
)

# ------------------------------------------------------------------------------
# Figure: GAM smooth of ideology (with political interest as control)
png("Figure_PI1.png", width = 800, height = 600)
plot(
  pi.gam_m1,
  select = 1,
  seWithMean = TRUE,
  shift = coef(pi.gam_m1)[1],
  shade = TRUE,
  xlab = "Ideology",
  ylab = "Populist Attitudes",
  cex.lab = 1.2,
  cex.axis = 1,
  rug = TRUE
)

dev.off()

# Figure: GAM smooth of extremism (with political interest as control)
png("Figure_PI2.png", width = 800, height = 600)
plot(
  pi.gam_m2,
  select = 1,
  seWithMean = TRUE,
  shift = coef(pi.gam_m2)[1],
  shade = TRUE,
  xlab = "Ideological Extremism",
  ylab = "Populist Attitudes",
  cex.lab = 1.2,
  cex.axis = 1,
  rug = TRUE
)

dev.off()

# ------------------------------------------------------------------------------
# Appendix E: Robustness Tests 
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# (1) – Data Preparation
df3 <- df %>%
  select(pop_add_m, pop_gz_mc, pop_mult_mc, pop_add,
         ideo, ext2,
         age, gen, edu, empl,
         corrupt_ind, econ_ind, pol_int,
         socecon_status,
         country_name) %>%
  mutate(
    country_name = factor(country_name),
    socecon_status = case_when(
      socecon_status == 1 ~ "White Collar",
      socecon_status == 2 ~ "Worker",
      socecon_status == 3 ~ "Farmer",
      socecon_status == 4 ~ "Self-employed",
      TRUE ~ NA_character_
    ),
    socecon_status = relevel(factor(socecon_status), ref = "White Collar")
  ) %>%
  filter(!(ideo == 5 & pol_int %in% c(3, 4)))  # Drop ambiguous center & uninterested

# ------------------------------------------------------------------------------
# (2) – Re-estimated Model 1 (Ideology): OLS Linear and Quadratic
# Filter Data: Exclude respondents with ideology == 5 and low interest (3, 4)
df2.2 <- df2 %>%
  filter(!(ideo == 5 & pol_int %in% c(3, 4)))
# Linear Model
re.ols_m1.3.lin <- lm(
  pop_add_m ~ ideo + age + gen + edu + empl +
    corrupt_ind + econ_ind + factor(country_name),
  data = df2.2
)

# Quadratic Model
re.ols_m1.3.quad <- lm(
  pop_add_m ~ ideo + I(ideo^2) + age + gen + edu + empl +
    corrupt_ind + econ_ind + factor(country_name),
  data = df2.2
)

# Clustered SEs by Country
re.ols_m1.3.lin_robust_clustered <- lmtest::coeftest(
  re.ols_m1.3.lin,
  vcov = sandwich::vcovCL,
  type = "HC1",
  cluster = ~df2.2$country_name
)

re.ols_m1.3.quad_robust_clustered <- lmtest::coeftest(
  re.ols_m1.3.quad,
  vcov = sandwich::vcovCL,
  type = "HC1",
  cluster = ~df2.2$country_name
)

# Table Output – Linear & Clustered
stargazer::stargazer(
  re.ols_m1.3.lin, re.ols_m1.3.lin_robust_clustered,
  type = "text",
  out = "Table_E2_2.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("ideo", "age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"
)

# Table Output – Quadratic & Clustered
stargazer::stargazer(
  re.ols_m1.3.quad, re.ols_m1.3.quad_robust_clustered,
  type = "text",
  out = "Table_E2_3.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("ideo", "age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"
)

# Plot: OLS Quadratic
p_re.ols_m1.3.quad <- ggeffects::ggpredict(
  re.ols_m1.3.quad, terms = "ideo",
  vcov.fun = "vcovCL", vcov.type = "HC1",
  vcov.args = list(cluster = df2.2$country_name)
)

png("Figure_E2c.png", width = 800, height = 600)
plot(p_re.ols_m1.3.quad) + geom_rug(sides = "b") +
  labs(x = "Ideology", y = "Populist Attitudes") +
  theme_bw() +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 13), axis.title = element_text(size = 16))
dev.off()

# ------------------------------------------------------------------------------
# (3) – Re-estimated Model 2 (Extremism): OLS Linear and Exponential
# Linear Model
re.ols_m2.3.lin <- lm(
  pop_add_m ~ ext2 + age + gen + edu + empl +
    corrupt_ind + econ_ind + factor(country_name),
  data = df2.2
)

# Exponential Model
re.ols_m2.3.exp <- lm(
  pop_add_m ~ ext2 + exp(ext2) + age + gen + edu + empl +
    corrupt_ind + econ_ind + factor(country_name),
  data = df2.2
)

# Clustered SEs by Country
re.ols_m2.3.lin_robust_clustered <- lmtest::coeftest(
  re.ols_m2.3.lin,
  vcov = sandwich::vcovCL,
  type = "HC1",
  cluster = ~df2.2$country_name
)

re.ols_m2.3.exp_robust_clustered <- lmtest::coeftest(
  re.ols_m2.3.exp,
  vcov = sandwich::vcovCL,
  type = "HC1",
  cluster = ~df2.2$country_name
)

# Table Output – Linear & Clustered
stargazer::stargazer(
  re.ols_m2.3.lin, re.ols_m2.3.lin_robust_clustered,
  type = "text",
  out = "Table_E2_5.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("ext2", "age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"
)

# Table Output – Exponential & Clustered
stargazer::stargazer(
  re.ols_m2.3.exp, re.ols_m2.3.exp_robust_clustered,
  type = "text",
  out = "Table_E2_6.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("ext2", "age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001"
)

# Plot: OLS Exponential
p_re.ols_m2.3.exp <- ggeffects::ggpredict(
  re.ols_m2.3.exp, terms = "ext2",
  vcov.fun = "vcovCL", vcov.type = "HC1",
  vcov.args = list(cluster = df2.2$country_name)
)

png("Figure_E2_OLS.png", width = 800, height = 600)
plot(p_re.ols_m2.3.exp) + geom_rug(sides = "b") +
  labs(x = "Ideological Extremism", y = "Populist Attitudes") +
  theme_bw() +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 13), axis.title = element_text(size = 16))
dev.off()

# ------------------------------------------------------------------------------
# E.3 – Populist Attitudes (Multiplicative)
# ------------------------------------------------------------------------------
# GAM: Ideology (Multiplicative PA)
re.gam_m1_pa.mc <- gam(
  pop_mult_mc ~ s(ideo) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df3
)

# GAM: Extremism (Multiplicative PA)
re.gam_m2_pa.mc <- gam(
  pop_mult_mc ~ s(ext2) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df3
)

# Diagnostics
gam.check(re.gam_m1_pa.mc)  # p-value for s(ideo) = 0.34
gam.check(re.gam_m2_pa.mc)  # p-value for s(ext2) = 0.09
concurvity(re.gam_m1_pa.mc, full = TRUE)  # Worst concurvity ~0.15
concurvity(re.gam_m2_pa.mc, full = TRUE)  # Worst concurvity ~0.16

# Model Summaries
summary(re.gam_m1_pa.mc)
summary(re.gam_m2_pa.mc)

# Table Output
stargazer::stargazer(
  re.gam_m1_pa.mc, re.gam_m2_pa.mc,
  type = "text",
  out = "Table_E3.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes (Multiplicative)",
  keep = c("age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001",
  notes.append = FALSE
)

# Visualization – SM Figure 5 & 6
png("Figure_E3_GAM.png", width = 800, height = 600)
plot(re.gam_m1_pa.mc, select = 1, seWithMean = TRUE, shift = coef(re.gam_m1_pa.mc)[1],
     shade = TRUE, xlab = "Ideology", ylab = "Populist Attitudes (Multiplicative)",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()

png("Figure_E3_GAM_ext.png", width = 800, height = 600)
plot(re.gam_m2_pa.mc, select = 1, seWithMean = TRUE, shift = coef(re.gam_m2_pa.mc)[1],
     shade = TRUE, xlab = "Ideological Extremism", ylab = "Populist Attitudes (Multiplicative)",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()

# ------------------------------------------------------------------------------
# E.4 – Populist Attitudes (Goertz Approach)
# ------------------------------------------------------------------------------
# GAM: Ideology (Goertz PA)
re.gam_m1_pa.gz <- gam(
  pop_gz_mc ~ s(ideo) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df3
)

# GAM: Extremism (Goertz PA)
re.gam_m2_pa.gz <- gam(
  pop_gz_mc ~ s(ext2) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df3
)

# Diagnostics
gam.check(re.gam_m1_pa.gz)  # p-value for s(ideo) = 0.79
gam.check(re.gam_m2_pa.gz)  # p-value for s(ext2) = 0.26

concurvity(re.gam_m1_pa.gz, full = TRUE)  # Worst concurvity ~0.15
concurvity(re.gam_m2_pa.gz, full = TRUE)  # Worst concurvity ~0.16

# Model Summaries
summary(re.gam_m1_pa.gz)
summary(re.gam_m2_pa.gz)

# Table Output
stargazer::stargazer(
  re.gam_m1_pa.gz, re.gam_m2_pa.gz,
  type = "text",
  out = "Table_E4.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes (Goertz)",
  keep = c("age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001",
  notes.append = FALSE
)

# Visualization – SM Figure 7 & 8
png("Figure_E4_GAM_ideo.png", width = 800, height = 600)
plot(re.gam_m1_pa.gz, select = 1, seWithMean = TRUE, shift = coef(re.gam_m1_pa.gz)[1],
     shade = TRUE, xlab = "Ideology", ylab = "Populist Attitudes (Goertz)",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()

png("Figure_E4_GAM_extremism.png", width = 800, height = 600)
plot(re.gam_m2_pa.gz, select = 1, seWithMean = TRUE, shift = coef(re.gam_m2_pa.gz)[1],
     shade = TRUE, xlab = "Ideological Extremism", ylab = "Populist Attitudes (Goertz)",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()

# ------------------------------------------------------------------------------
# E.5 – Populist Attitudes (without Manichaeism)
# ------------------------------------------------------------------------------
# GAM: Ideology (w/o MC)
re.gam_m1_pa.woMC <- gam(
  pop_add ~ s(ideo) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df3
)

# GAM: Extremism (w/o MC)
re.gam_m2_pa.woMC <- gam(
  pop_add ~ s(ext2) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    factor(country_name),
  method = "REML",
  data = df3
)

# Diagnostics
gam.check(re.gam_m1_pa.woMC)  # p = 0.48
gam.check(re.gam_m2_pa.woMC)  # p = 0.80
concurvity(re.gam_m1_pa.woMC, full = TRUE)  # max ~ 0.14
concurvity(re.gam_m2_pa.woMC, full = TRUE)  # max ~ 0.16

# Summaries
summary(re.gam_m1_pa.woMC)
summary(re.gam_m2_pa.woMC)

# Table Output
stargazer::stargazer(
  re.gam_m1_pa.woMC, re.gam_m2_pa.woMC,
  type = "text",
  out = "Table_E5.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes (w/o MC)",
  keep = c("age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001",
  notes.append = FALSE
)

# Visualization: SM Plot 09 & 10
png("Figure_E5_GAM_ide.png", width = 800, height = 600)
plot(re.gam_m1_pa.woMC, select = 1, seWithMean = TRUE, shift = coef(re.gam_m1_pa.woMC)[1],
     shade = TRUE, xlab = "Ideology", ylab = "Populist Attitudes (w/o MC)",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()

png("Figure_E5_GAM_ext.png", width = 800, height = 600)
plot(re.gam_m2_pa.woMC, select = 1, seWithMean = TRUE, shift = coef(re.gam_m2_pa.woMC)[1],
     shade = TRUE, xlab = "Ideological Extremism", ylab = "Populist Attitudes (w/o MC)",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()

# ------------------------------------------------------------------------------
# E.6 – Populist Attitudes with Socioeconomic Status as Control
# ------------------------------------------------------------------------------
# GAM: Ideology + Socioeconomic Status
re.gam_m1_pa.socecon <- gam(
  pop_add_m ~ s(ideo) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    socecon_status +
    factor(country_name),
  method = "REML",
  data = df3
)

# GAM: Extremism + Socioeconomic Status
re.gam_m2_pa.socecon <- gam(
  pop_add_m ~ s(ext2) +
    age + gen + edu + empl +
    corrupt_ind + econ_ind +
    socecon_status +
    factor(country_name),
  method = "REML",
  data = df3
)

# Diagnostics
gam.check(re.gam_m1_pa.socecon)  # p = 0.9
gam.check(re.gam_m2_pa.socecon)  # p = 0.3
concurvity(re.gam_m1_pa.socecon, full = TRUE)  # max ~ 0.13
concurvity(re.gam_m2_pa.socecon, full = TRUE)  # max ~ 0.12

# Summaries
summary(re.gam_m1_pa.socecon)
summary(re.gam_m2_pa.socecon)

# Table Output
stargazer::stargazer(
  re.gam_m1_pa.socecon, re.gam_m2_pa.socecon,
  type = "text",
  out = "Table_E6.txt",
  odd.ratio = TRUE,
  dep.var.labels = "Populist Attitudes",
  keep = c("age", "gen", "edu", "empl", "corrupt_ind", "econ_ind", "socecon_status", "Constant"),
  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
  star.char = c("+", "*", "**", "***"),
  ci = TRUE,
  align = TRUE,
  notes = "Odds with 95% Confidence Intervals in parentheses. + p<0.1; * p<0.05; ** p<0.01; *** p<0.001",
  notes.append = FALSE
)

# Visualization: SM Plot 11 & 12
png("Figure_E6_GAM_ideo.png", width = 800, height = 600)
plot(re.gam_m1_pa.socecon, select = 1, seWithMean = TRUE, shift = coef(re.gam_m1_pa.socecon)[1],
     shade = TRUE, xlab = "Ideology", ylab = "Populist Attitudes",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()

png("Figure_E6_GAM_ext.png", width = 800, height = 600)
plot(re.gam_m2_pa.socecon, select = 1, seWithMean = TRUE, shift = coef(re.gam_m2_pa.socecon)[1],
     shade = TRUE, xlab = "Ideological Extremism", ylab = "Populist Attitudes",
     cex.lab = 1.2, cex.axis = 1, rug = TRUE)
dev.off()


# Stop logging------------------------------------------------------------------
# Log session information
sessionInfo()
sink()